Thermal depinning of fluxons in discrete Josephson rings 
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We study the thermal depinning of single fluxons in rings made of Josephson junctions. Due to 
thermal fluctuations a fluxon can be excited from its energy minima and move through the array, 
causing a voltage across each junction. We find that for the initial depinning, the fluxon behaves as 
a single particle and follows a Kramers-type escape law. However, under some conditions this single 
particle description breaks down. At low values of the discreteness parameter and low values of the 
damping, the depinning rate is larger than the single particle result would suggest. In addition, for 
some values of the parameters the fluxon can undergo low-voltage diffusion before switching to the 
high- voltage whirling mode. This type of diffusion is similar to phase diffusion in a single junction, 
but occurs without frequency-dependent damping. We study the switching to the whirling state as 
well. 
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I. INTRODUCTION 

In past years many works have been devoted to the 
study of extended discrete nonlinear systems. On the 
one hand, it is important to deepen our knowledge of 
general properties of such systems since they often have 
application to many different physical situations. On the 
other hand, many physical systems are well described by 
nonlinear discrete models. In this field, the emergence 
of the concept of soliton for instance (in the continu- 
ous and its discrete counterparts) was paradigmatic. 1,2 
A well known example of model system supporting this 
type of nonlinear excitations is the discrete Sine-Gordon 
equation (also called Frenkel-Kontorova model) 

A Josephson-junction (JJ) array is by construction a 
discrete system made of interacting nonlinear solid state 
devices. J J arrays are an example of physical systems 
with great fundamental and technological interest which 
are well described by discrete nonlinear models From 
the experimental point of view, Josephson junctions are 
a privileged place to study solitons and to explore their 
possible applications^ 

Of the many different geometries for a JJ array the so- 
called Josephson ring (a set of JJ connected in parallel 
and closed forming a ring, see Fig. [T]) is well described by 
a discrete sine-Gordon equation and supports nonlinear 
discrete solitons or kinks, usually called fluxons in this 
context^ Thus, many studies of the discrete sine-Gordon 
equation or the role of kinks in nonlinear arrays have 
direct application when studying JJ rings. Conversely 
the study and modelization of J J rings allows exploration 
of new issues concerning the behavior of these nonlinear 
phenomena. 

Our model system is the so-called Josephson ring, a 
collection of JJ coupled in parallel and forming a ring 
(see figure [T]). When cooled below the superconducting 
critical temperature an integer number of magnetic flux 
quanta, fluxons, can be trapped in the ring. Then, the 
physical properties of the array are dominated by the 



presence of the fluxons in the system. The I-V curve of 
the array shows the mean voltage across the array as a 
function of a constant external current applied to every 
junction of the array. When current is applied, the array 
remains superconducting up to a certain critical value 
which defines the critical current of the array. In the 
presence of fluxons this current corresponds to the fluxon 
depinning current. At this current, the fluxon starts to 
move around the array and finite voltage is measured. 

In many cases the energy exchange between the system 
and the environment is relevant, so thermal fluctuations 
have to be considered^ Because of this, at finite temper- 
ature the value of the depinning current and the shape of 
the I-V curve can be strongly affected by thermal fluctu- 
ations. Motivated by recent experiments on the thermal 
depinning and dynamics of fluxons (kinks) in small rings 
made of 9 junctions^ in this work we numerically ex- 
plore some of these issues. In addition, in some cases 
the fluxons can be understood as particles on a substrate 
periodic potential. 

The main object of this paper is to study numerically 
the thermal depinning of a single fluxon in a JJ ring 
and compare it with numerical simulations and analyti- 
cal predictions for the case of a single particle. We have 




FIG. f: Left: Scheme of the JJ ring. Lines are for super- 
conducting wires and crosses for Josephson junctions. Right: 
Phase configuration profile (top) and magnetic flux (bottom) 
for a fluxon in J J ring. 
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found excellent agreement in many cases. However, un- 
der some conditions the single particle description fails. 
In addition, for some values of the parameters the fluxon 
can undergo low-voltage diffusion before switching to the 
high- voltage whirling mode. This type of diffusion is sim- 
ilar to phase diffusion in a single underdamped junction, 
but occurs without frequency-dependent damping. 



II. EQUATIONS 

Josephson junctions are made of two superconduct- 
ing materials separated by a thin insulating barrier. 
Driven by an external current, this system behaves as 
a solid-state nonlinear oscillator and is modeled by the 
same dynamical equations that describe a driven pen- 
dulum: 7 i = (p + r^J + sin<y9 + £(t). Here ip, the vari- 
able that describes the behavior of the junction, is the 
gauge-invariant phase difference of the superconducting 
order parameter at both sides of the junction. In this 
equation current is normalized by the junction critical 
current I c and time by the junction plasma frequency 
ujp = ^/2nI c /QoC ($0 = h/2e is the magnetic flux quan- 
tum and C the junction capacitance). Y is an impor- 
tant parameter which measures the dissipation in the 
system (Y — ^/<5>o/27r/ c Ci? 2 , with R the effective re- 
sistance of the junction). The last term, £(r), describes 
the effect of thermal noise in the dynamics (Johnson cur- 
rent noise) and satisfies (£(t)) = and (£(t)£(t')) = 
2YT5(t — t') where we use T for a normalized temper- 
ature T = fcsToxp/ Ej (with Ej the Josephson energy 
Ej = ^qIc/2tt). The normalized dc voltage v which 
gives the response of the system to the external cur- 
rent is defined by v — Vde/hR = {§o/2nI c R){dp/dt) = 
Y(d(p/dr}. 

As previously stated, the JJ ring consists of a series of 
individual junctions connected in parallel. Such system 
can be though of as a series of coupled pendula. Following 
the usual model for the system, the equations for an array 
made of N coupled junctions driven the same external 
current are given by£ 

ipj + Yipj + sin ipj + £j (r) = X((pj+i —2<pj + (pj_i)+i (1) 

Index j denotes different junctions and run from I to N. 
The new parameter A accounts for the coupling between 
the junctions which, in the framework of this model, oc- 
curs between neighbors and has an inductive character 
A = $o/2ttI c L, with L the self inductance of every cell 
in the array. Boundary conditions are defined by the 
topology of the array (here we consider circular arrays) 
and the number M of trapped fluxons in the system: 
fj+N = <Pj + 2nM. 

In this article we will consider the case of a single 
fluxon. We have studied different sizes for the array, but 
here we will present results for an array made of 9 junc- 
tions, similar to those being experimentally studied. We 
will also study different values of A and restrict our in- 



terest to underdamped arrays biased by a dc current in 
a broad range of temperatures. 



III. RESULTS 

In this section we are going to present numerical simu- 
lations of the dynamics of one fluxon in a Josephson ring 
and one particle in a periodic potential. We will also 
show numerical calculations from single particle thermal 
escape theory. 



A. I-V curves: damping regimes 

Figure [2] shows single I-V curves for one fluxon in a 9 
junctions Josephson ring with A = 0.4. In this case the 
width of the fluxon is close to 2, so it is well localized in 
the array and discreteness effects are important*^ Curves 
were simulated at three different values of damping and 
five temperatures. Current was increased from zero to 
some maximum at an average ramp equal to | x I0~ 7 in 
normalized units. 

Let us look first at the T = curves. If we start at 
zero bias, in all the cases the fluxon is pinned to the array 
up to the so-called depinning current i^ ep is reached (for 
A = 0.4, ij cp ~ 0.155). Above this current very different 
I-V characteristics are observed depending on the value 
of the damping. 

At small damping the system switches from the v = 
state to a high-voltage ohmic state (v = i) where all the 
junctions rotate uniformly (whirling branch). The damp- 
ing is so small that when the fluxon moves through the 
array it excites all the junctions to the high- voltage state. 
In this voltage state the fluxon is totally delocalized in 
the array. For clarity, we have shown only curves for 
increasing current. If current is decreased from the high- 
voltage state the junctions retrap at a small value of the 
current (mostly defined by Y). The curve is hysteretic 
and shows bistability for a wide range of currents. 

At intermediate values of damping the dynamics is 
much more complex. Now the I-V curve shows a low- 
voltage region dominated by a series of steps which cor- 
respond to resonances between the fluxon velocity and 
the linear modes of the array. These resonances has been 
the object of great attention in the pasti£ ' 12 i 13 ' 14 i 15 ' 16 i 17 
For these currents the fluxon moves around the ring in 
a localized manner. This regime persists up to a given 
value of the current for which the fluxon reaches a high 
velocity and all the junctions switch to the high-voltage 
part of the curve. At moderate damping the IV curve can 
also be multistable with hysteresis loops on every step. 

At high values of the damping dissipation governs the 
dynamics, multistability disappears and the voltage in- 
creases from zero without discontinuities and jumps as 
soon as current reaches the depinning value. In this part 
of the curve a localized fluxon is moving around the ring 



3 



T=0 

T^O.002 
T=0.01 
T=0.02 
T=0.1 




0,001 0,002 

voltage 



0,02 0,04 

voltage 



FIG. 2: (Color online) I — V curves for one fluxon in a 9 junctions ring with A = 0.4 at T = 0.01 (left), T = 0.1 (middle) and 
r = 1.0 (right). Each figure shows 5 different temperatures (T = 0, 0.002, 0.01, 0.02, 0.1). 



and voltage is related with the fluxon velocity. For cur- 
rents close to 1 (in normalized units) it starts the tran- 
sient to another regime where the fluxon delocalizes and 
all the junctions rotate and contribute to the overall volt- 
age in the array£ 

Figure [2] also shows the dynamics of the array in the 
presence of thermal noise. As can be seen in the figure, 
the first thermal effect is that the fluxon depins at smaller 
currents. For small damping we find that if tempera- 
ture is high enough, noise also induces a fluxon diffusion 
branch, which we discuss later in the paper. At moder- 
ate damping, the low-voltage resonances are rounded and 
at high enough temperature voltage increases smoothly 
from zero to some value on the fluxon diffusion branch, 
and then switches to the high-voltage region of the 1-V 
curve. At high damping temperature causes a rounding 
of the curve. 

We will study how temperature affects the FV curve 
at low damping since this is the case for the experimen- 
tal system we are trying to model. Usually the depinning 
current is experimentally defined as the current for which 
measured voltage is above a certain threshold. We have 
followed the same definition in our simulations. This is a 
good definition in the low damping and low temperature 
regime where the system switches between two very dif- 
ferent voltage values so a threshold independent current 
is expected. However, the election of the threshold volt- 
age is not trivial. A small threshold can give problems 
at large temperatures where voltage fluctuations are also 
large. A large threshold is not a good choice since ignore 
possible low-voltage states like the fluxon diffusion one. 
If low-voltage states are present we distinguish between 



the depinning current idep and the switching current i sw , 
where depinning marks the end of the superconducting 
state and switching the transition to the high-voltage 
branch. At very high temperatures noise can reach the 
threshold level and first switching is suppressed. To study 
this issue we have used in our simulation three or five dif- 
ferent thresholds and compare results for all of them. 



B. (i dcp (T)) 

The depinnig current is a stochastic variable with a 
given probability distribution. We present results for the 
mean value of the depinning current (idep) and its stan- 
dard deviation a. Results were obtained after the numer- 
ical simulation of the dynamical equations of the system 
for 1000 samples. We used different values of damping 
r (typically from 0.001 to 0.1) and coupling A (usually 
from 0.2 to 1.0). The number of junctions in the array is 
N = 9. Another important parameter of the simulations 
is the current ramp; in our case this ramp change from 
one simulation to other but it is of the order of 10~ 7 . 

Figure [3] shows results for T — 0.01 and five different 
values of the coupling (A = 0.2, 0.4, 0.6, 0.8 and 1.0). The 
main physical properties of the fluxon in the array change 
importantly with the value of A. Thus, the zero tempera- 
ture depinning current of the array i° ep decreases a factor 
of 30 from A = 0.2 to A = 1.0 (see table below). In or- 
der to compare the five curves, in figure[3]both axis have 
been scaled by the value of the zero temperature depin- 
ning current for every case J£ We see that once scaled all 
curves are similar showing that in this range of parame- 
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FIG. 3: (Color online) Numerical calculation of normalized 
(idep) an d cr versus normalized T at V — 0.01 and 5 different 
values of A (A = 0.2, 0.4, 0.6, 0.8, 1.0). 



ters these results can be understood in a unified manner. 

In the figure we see that («dep) decreases to zero as 
effect of temperature. All the curves follow the same 
behavior but at high temperatures the small A curves 
slightly deviates from the others. For temperatures of 
the order of the barrier thermal fluctuations dominate 
the dynamics and the depinning current goes to zero. In 
fact, at high temperatures there is no a good definition 
of depinning current since different thresholds may give 
different results. With respect to the standard devia- 
tion, we can see that it grows with T 2 / 3 as predicted by 
standard thermal activation theory and reaches a maxi- 
mum at high temperatures, when (idep) has an inflection 
point, close to (idep) —> 0. Then the standard deviation 
decreases since all escape events happens in a narrow 
range of current values. 



C. The fluxon as a single particle. 

When studying the dynamics of one fluxon in the ar- 
ray it is very common to use the picture of this extended 
and collective object as a single particlei^ i 19 ' 20 This ap- 
proach has been extensively used in the past and, as we 
will see, it is very useful although not exact. 



TABLE I: Fluxon parameters at different values of A 
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£pn 


2 

LO PN 


in 




0.2 


0.77842 


0.77405 


0.5028 


0.38482 


0.4 


0.30974 


0.44775 


0.3459 


0.15435 


0.6 


0.12744 


0.23151 


0.2757 


0.06367 


0.8 


0.05539 


0.11721 


0.2363 


0.02767 


1.0 


0.02550 


0.06041 


0.2110 


0.012725 



Let us consider a new variable representing the center 
of masses of the fluxon or the position of the fluxon in 
the array. Then, in the simplest approach, the dynamics 
of a fluxon in a ring can be approach by the dynamics 
of a driven, damped, massive particle experiencing a si- 
nusoidal substrate potential (Peierls-Nabarro potential) 
and subjected to thermal fluctuations: 

mX + TmX + ^ op sin X = i + £(t) (2) 

where 

(£(t)) = and (f (t)£(t')> = 2mTT5{r - r') (3) 

In this simple approach we are neglecting for instance 
the spatial dependence of the mass, effective damping due 
to the other degrees of freedom of the system and higher 
order terms in the expansion of the substrate potential 
for the fluxon. 

Table U gives a relation of numerically computed values 
of some of the parameters of the fluxon and its effective 
potential in the single particle picture: £pn i s the zero 
current potential barrier; Up N is the squared frequency 
for small amplitude oscillations of the fluxon around equi- 
librium; m is the fluxon effective mass at rest (computed 
as in — _Epn/2o->p N ) and i^ ep the depinning current. For a 
perfect sinusoidal potential we should get i^ cp = Ep^/2. 
The exact results are close to it. 

Figures [4] and [5] show for A = 0.4 and A = 0.8 respec- 
tively (in both cases T — 0.01) the comparison between 
the results for the fluxon in the array, numerical simula- 
tions of the depinning of a single particle in a sinusoidal 
potential (Eq. [2]) and a theoretical calculation based on 
analytical results for the thermal activation rate of par- 
ticles in sinusoidal potentials in the low damping regi- 
men ,2ii22i23i24i25i26 ^ ^ figure we plot the result com- 
puted from the Buttiker, Harris and Landauer equation 
for escape rate at a = 1 [Eq. (3.11) in reference^ 3 .]. We 
have checked that a very close result (indistinguishable 
at the scale of the figure) is got when using a = 1.45 — 
or the Melnikov and Meshov theory^ 

We can see that the single particle simulations repro- 
duce the results for the fluxon in an excellent way at 
this value of the damping in all the temperature range 
although a small deviation in a range of temperatures is 
observed for A = 0.4. Theoretical estimations disagree 
at high values of T since escape rate equations were ob- 
tained in the infinite barrier limit (Eb 3> fcsT). 
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FIG. 4: (i dcp ) and a versus T at V = 0.01 and A = 0.4 for 
the fluxon and a single particle in a periodic potential and 
comparison with the theoretical prediction. 
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FIG. 5: (idcp) and a versus T at Y = 0.01 and A = 0.8 for 
the fluxon and a single particle in a periodic potential and 
comparison with the theoretical prediction. 



The agreement shown in Figures |4] and [5] does not oc- 
cur for other values of coupling and damping. For in- 
stance, for A = 0.4 and T = 0.001 (figure [9] below) a 
deviation of the simulated curve with respect to the the- 
oretical prediction is found. To study further such result 
we have done numerical simulations at fixed T for differ- 
ent values of T and for A = 0.4 (T = 0.01) and A = 0.8 
(T = 0.0018) ,22 Results are shown in figures © and [7] 
There we can see that for A = 0.4 the fluxon results devi- 
ate importantly from the expected for the single particle 
(or the theory) when damping is decreased. However this 
is not the case for A = 0.8. In this respect it seems to be 
important the degree of discreteness of the system. Such 
degree is measured by the coupling parameter A (high A 
approach the continuum limit of the system, and small 
A increases the discreteness effects). At small values of 
A effects due to other degrees of freedom are more im- 
portant and if damping is small such excitations persist 
longer in the system. 



D. Fluxon diffusion 

In figure [5] we have seen that at low damping and high 
enough temperature a low-voltage branch appears in the 
/ — V curves before the escape to the full running state. 



In this low-voltage state, the transport of the fluxon oc- 
curs through a series of noise-induce 2ir phase slips (or 
2?t/N depending on the definition of the fluxon center of 
masses) where every jump corresponds to a fluxon which 
advances one cell in the array. Such state can not be 
understood in terms of the single particle picture. In 
analogy with the phase diffusion that occurs in a single 
junction, we label this mode of transport fluxon diffusion. 
In spite of the fact of we are in the low damping regime, 
the fluxon is able to travel along the ring without excit- 
ing the whirling branch. Remarkably this is a thermally 
excited state and it is not seen at low temperatures. 

In figure [8] we show the time evolution of the phase 
of one junction (junction 1) and a phase associated with 
the center of mass of the fluxon defined as tp — y\ cpj 
to allow a better comparison. Current values have been 
chosen in the low-voltage branch The random jumps of 
2tt for the junction or 2-k/N for the fluxon phase ip are 
easily observed. 

We have also done numerical experiments to simulate 
the jump from the diffusion branch to the full whirling 
state. These simulations are done with a higher voltage 
threshold, but are otherwise similar to the previous ex- 
periments. Results for T — 0.01 and 0.001 are shown 
in figure O here A = 0.4. We see that theory and sin- 
gle particle results agree quite well in all the range. At 
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FIG. 6: (idop) as a function of r for A = 0.4 at T = 0.01. 



0,025 



0,020 - 



0,015 



1E-5 



1E-4 



1E-3 



0,01 



FIG. 7: (i dc p) as a function of T for A = 0.8 at T = 0.0018. 



high temperature thermal fluctuations dominate, are of 
the order of the barrier, and theoretical results do not 
apply. We also see that the standard deviation increases 
with temperature following the expected law up to a cer- 
tain value where reaches a maximum and then decreases 
when (idop) is close to zero. 

However, in figure [9] we can also see that for the fluxon 
curve there exists a value of T such that at higher tem- 
peratures, the value of the current for which the array 
switches to the whirling branch is temperature indepen- 
dent. This temperature shows the emergence of a fluxon 
diffusion branch in the I-V curve. Looking at I— V curves 
it appears that the switching current value is defined by 
some value of the fluxon velocity. Comparing the stan- 
dard deviation curve in figure [9] to the similar one with 
the lower threshold in figure |U as expected, we can see 
that a peak occurs much earlier in temperature for the 
jump from the diffusion state. 
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FIG. 8: A = 0.4, T = 0.01 and T = 0.1. Top: I-V curve 
showing the small voltage fluxon diffusion branch. Medium 
and bottom: time evolution of the phase of junction 1 (black 
line) and the phase of the center of masses of the fluxon (grey 
line) divided by 2tt at i = 0.005, 0.01, 0.02, 0.05 and 0.07. 



IV. DISCUSSION AND CONCLUSIONS 

We have studied the thermal depinning of fluxons in 
small Josephson rings at small values of damping. At zero 
temperature as current is increased the system switches 
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FIG. 9: (Color online) (w) and a as a function of T at 
A = 0.4, for two different values of the damping V — 0.01 
and r = 0.001. We show plots for the theoretical prediction 
(solid line), the particle (open symbols) and the fluxon (solid 
symbols) . 



from a superconducting zero voltage state to a resistive 
state where v = i. This happens at the so called de- 
pinning current. Beyond this current there is not static 
configuration for one fluxon in the array, and the fluxon 
starts to move. Due to the low value of the damping when 
fluxon goes through the array causes all the junctions to 
switch to a high-voltage state. Then all the junctions 
do the same but with a phase difference that accounts 
for the presence of one quantum of flux homogeneously 
distributed along the whole array. 

Due to thermal fluctuations, in an experiment the 
value of the measured depinning current changes from 
one I-V to another and only a probability distribution 
function has sense. This function is usually character- 
ized by its mean value and its standard deviation. The 
main object of this paper has been to numerically study 
how these observables behave for different system param- 
eters (coupling A, damping T and temperature T)/2S We 
also have compared these results with numerical simula- 
tions and theoretical estimations for the depinning of a 
single particle in a sinusoidal potential. 

As expected, the mean value of the depinning or the 
switching current decreases as temperature is raised. 
At low temperatures the standard deviation follows the 



usual T 2 / 3 law. At higher temperatures the cr(T) func- 
tion reaches a peak, that we can identify with the points 
on the (idep(T)) or (? SW (T)) curves at which the curvature 
changes sign (inflection point). 

Roughly speaking our results show that the depinning 
of the fluxon can be understood in terms of the stochastic 
dynamics of a single particle in a tilted sinusoidal poten- 
tial. However, we have seen some unexpected effects that 
we attribute to discreteness. Then for the case of small 
coupling (A = 0.4 and smaller) we have seen an increasing 
deviation of the fluxon depinning behavior from our ex- 
pectations from the single particle picture. Arrays with 
a larger coupling are closer to the continuous limit and 
discreteness effects are smaller; here the single particle 
picture works much better. 

The other effect we have observed is the emergence at 
small damping of a low-voltage thermally excited state 
that we call fluxon diffusion. In such cases the zero tem- 
perature I-V curve does not show any resonance or low- 
voltage state and the system switches from zero voltage 
to the high-voltage branch. However, after some tem- 
perature a low-voltage branch is observed. The system 
first switches from zero to this branch and then to the 
high-voltage state. At higher temperatures the system 
first continuously increases voltage from zero (then is not 
clear how to define idep) and at larger currents switches 
from the low-voltage state to the high-voltage one. We 
have also seen that the values of damping and current for 
which this behavior is observed depends importantly in 
A and also in the number of junctions in the array, N£Q 

An important point of our work has been to compare 
our numerical results with results based on the single 
particle picture. The main conclusion is that for most 
of the cases this picture gives a good estimation of the 
fluxon depinning current. In fact, in an experimental 
case, where the different parameters (mainly A, T and 
I c ) are known with some imprecision will be difficult to 
identify deviations from the expected behavior. In ad- 
dition, there are some points difficult to address: the 
effective one-dimensional potential for the fluxon in the 
array (Peierls-Nabarro potential) is not purely sinusoidal, 
and the value of the fluxon mass is not constant since it 
depends on the fluxon position and the current value. 
We are also neglecting all the system degrees of freedom 
except one and we know that in some cases this is not 
valid: for instance, to understand resonant steps which 
are due to coupling between the fluxon velocity and the 
linear waves of the discrete array or to understand the 
fluxon diffusion branch. We have also considered the ex- 
pression for the escape rate in a multidimensional case. 
In this expression usually the attempt frequency depends 
on the frequency of all the stable modes in the minimum 
and the saddle. ~» We have computed such numbers and 
check that the maximum error is smaller of 7% (and oc- 
curs for A = 0.125). 

Our numerical results for the single particle agree 
pretty well the predictions from Kramers theory for es- 
cape rate except for some limits. Disagreement at high 
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temperatures is expected since theoretical expressions 
are computed in the infinite barrier limit of the system 
(Eb ~S> kT). This limit is not fulfilled at high temper- 
atures. This is also true at small temperatures, where 
most of the escape events occur at currents very close to 
*dep where the barrier is also very small. We have checked 
that the E^/kT ratio in this case is also small. To fin- 
ish we have to mention the unexpected disagreement at 
small values of T. We are currently study further such 
results. 

In the single-particle picture, or the RCSJ model for 
a single junction, it has become generally accepted that 
diffusion cannot coexist with hysteresis. This was elu- 
cidated nicely in Kautz and Martinis^ through phase 
space arguments. Simply put, if the value of applied cur- 
rent is sufficient to allow a stable running state to coexist 
with the fixed points of zero voltage, the basins of at- 
traction for that running state necessarily separates the 
basins of attraction for any two neighboring fixed points. 
Phase jumps between two fixed points are thus forbid- 
den, as the system must first pass through the basin of 
attraction for the running state. While we have shown 
in previous sections that the initial escape of the fiuxon 
from its minima can be explained by thermal activation 
of a single particle, the fiuxon diffusion state in the I-V 
curves of figure [2] cannot be explained in a similar way. 

In the original single-junction experiments on phase 
diffusion, the coexistence of phase diffusion and hystere- 
sis was explained by the presence of frequency-dependent 
damping from the junction leads. A simple model of 
frequency-dependent damping is a series-RC circuit in 
parallel with the junction, which adds an extra degree of 
freedom to the phase space for the junction dynamics. 
This extra dimension resolves the above-mentioned issue 
regarding the non-overlapping basins of attraction. A 
main result of the paper is the observation of fiuxon dif- 
fusion in our simulations without frequency dependent- 



damping, which has not been included in our equations. 
Instead of the extra dimension introduced by frequency- 
dependent damping, fiuxon diffusion must occur because 
of the additional degrees of freedom from the multiple 
junctions in the array. We plan to explore the physics of 
this new diffusion mechanism in future experiments and 
simulations. 

When studying the behavior of the system at differ- 
ent temperatures we have seen that the mean value of 
the distribution of the switching current from the fiuxon 
diffusion branch is almost constant and the standard de- 
viation is very small. Comparing the standard deviation 
curve in figure [5] to the similar one with the lower thresh- 
old in figure [H we also see that a peak occurs much ear- 
lier in temperature for the jump from the diffusion state. 
This peak in the standard deviation is also reminiscent of 
experiments on single junctions ; 32 ' 33 ^ 4 where a peak in 
the standard deviation indicated the collapse of thermal 
activation and the onset of diffusion. 

To finish we want to mention that to our knowledge 
there are not experimental or numerical results studying 
systematically the thermal escape of fluxons or solitons 
in discrete arrays. However we want here to mention 
the work by A. Wallraff et al on vortices in long JJM- 
With respect with theoretical advances we want also to 
cite recent work 3 ^ where major differences between the 
macroscopic quantum tunnelling in JJ from tunnelling of 
a quantum particle are reported. 
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